#$maindir="/data1/bsi/BORA_processing/devel/data_HTHGU/illumina1/MayoGAP/mayogap_emerge/raw/clean/beagle/run_beagle/chr22/1000samples";
#$maindir="/data1/bsi/BORA_processing/devel/data_HTHGU/illumina1/MayoGAP/mayogap_emerge/raw/clean/beagle/run_beagle/chr22";
$maindir="/data1/bsi/BORA_processing/devel/data_HTHGU/illumina1/MayoGAP/mayogap_emerge/raw/clean/beagle/run_beagle/chr22";
$markerfile = "$maindir/whole_default/markers.txt";
#$markerfile = "/data1/bsi/BORA_processing/devel/data_HTHGU/illumina1/MayoGAP/mayogap_emerge/raw/clean/beagle/run_beagle/chr22/chr22_map.txt";
#$r2file = "$maindir/whole_default/beagle_out_chr22.geno.bgl.gz.r2";
$r2file = "$maindir/whole_default/beagle_out_chr22.geno.bgl.gz.r2";
$num_subdirec=4;
#first file shud be default run
#middle file shud be the edges
#last file is the random seed run
@files = ("$maindir/whole_default/beagle_out_chr22.geno.bgl.gz.gprobs","$maindir/quaterquatermb/whole_gprobs","$maindir/halfquatermb/whole_gprobs","$maindir/quatermb/whole_gprobs","$maindir/halfmb/whole_gprobs","$maindir/1mb/whole_gprobs","$maindir/2mb/whole_gprobs","$maindir/4mb/whole_gprobs","$maindir/whole_random2/beagle_out_chr22.geno.bgl.gz.gprobs");
#creating the file names
#$line="";
@files_edge = ("quaterquatermb","halfquatermb","quatermb","halfmb","1mb","2mb","4mb","seed");
@filenames = ("sumdiff","meandiff","meddiff","per1diff","per2diff","per5diff");
@seg_edges=("quaterquatermb","halfquatermb","quatermb","halfmb","1mb","2mb","4mb");
#creating the map file
for($i=1;$i<@files;$i++)
{
	print "dealing with $files[$i]\n";
	$line="";
	#creating the header
	for($j=0;$j<@filenames;$j++)
        {
                $temp = $filenames[$j].'.'.$files_edge[$i-1];
                $line = $line.$temp."\t";
        }
	$line =~ s/\t$//g;
	$filemap = "$maindir/$files_edge[$i-1]".'_temp.txt';
	print "FILEMAP $filemap\n";
	open(MAP,">$filemap") or die "not able to write $filemap\n";
	print MAP $line."\n";
	open(REF,$files[0]) or die "no file exists $files[0] $!\n";
	open(ALT,$files[$i]) or die "no file exists $files[$i] $!\n";
	$ref = <REF>;
	$alt = <ALT>;
	#if headers of dosage file is not equal
	if($ref ne $alt)
	{
		die "samples are not equal or not in the same order $files[0] & $files[$i]\n";
	}
	while($ref = <REF>)
	{
		$alt = <ALT>;
		chomp($ref);
		chomp($alt);
		@ref = split(" ",$ref);
		@alt = split(" ",$alt);
		if($ref[0] ne $alt[0])
		{
			die " snps not equal $ref[0] $alt[0] $files[0] $files[$i]\n";
		}
		@dosagediff=();
		for($j=3;$j<@ref;$j++)
		{
			$first = $j;
			$j++;
			$second = $j;
			$j++;
			$val = abs((2*$ref[$first]+$ref[$second])-(2*$alt[$first]+$alt[$second]));
			push(@dosagediff,$val);
		}
		@values = ();
		push(@values,sum(\@dosagediff));
		push(@values,mean(\@dosagediff));
		push(@values,median(\@dosagediff));
		push(@values,per1diff(\@dosagediff));
		push(@values,per2diff(\@dosagediff));
		push(@values,per5diff(\@dosagediff));
		$line = join("\t",@values);
		print MAP $line."\n";
	}
}
#combining all the files
#open(BUFF,$markerfile) or die "no file found $markerfile\n";
#open(BUFF0,$r2file) or die "no file found $r2file\n"; 
system("paste $markerfile $r2file > $maindir/map_temp.txt");
$fi = 'awk \'{if($1==$5)print $1"\t"$2"\t"$3"\t"$4"\t"$6}\'  '."$maindir/map_temp.txt> $maindir/map_temp1.txt";
system($fi);
system("mv $maindir/map_temp1.txt $maindir/map_temp.txt");
open(BUFF0,"$maindir/map_temp.txt") or die "no file found $maindir/map_temp.txt\n";

open(WRBUFF,">$maindir/map.txt") or die "not able to write\n";
#open(WRBUFF1,">$maindir/segdat.txt") or die "not able to write\n";
#merging the map header
print WRBUFF "rsid\tbp\talleleA\talleleB\tr2";
for($i=1;$i<@files;$i++)
{
        print "dealing with $files[$i]\n";
        $filemap = "$maindir/$files_edge[$i-1]".'_temp.txt';
	$buff="BUFF$i";
	open($buff,$filemap) or die "no file exists $filemap\n";
	$line=<$buff>;
	chomp($line);
	print WRBUFF "\t$line";
}
print WRBUFF "\n";
#merging the files
while($line=<BUFF0>)
{
	chomp($line);
	print WRBUFF $line;
	for($i=1;$i<@files;$i++)
	{
		$buff="BUFF$i";
		$line=<$buff>;
		chomp($line);
		print WRBUFF "\t$line";
	}
	print WRBUFF "\n";
}
#removing the temp files
system("rm $maindir/map_temp.txt ");
for($i=1;$i<@files;$i++)
{
	$filemap = "$maindir/$files_edge[$i-1]".'_temp.txt';
	system("rm $filemap");
}
#creating the seg file
open(WRBUFF,">$maindir/segdat.txt") or die "not able to write segdat.txt\n";
print WRBUFF "break.";
for($i=0;$i<@seg_edges;$i++)
{
	print WRBUFF "\tseg".'.'."$seg_edges[$i]"
}
print WRBUFF "\n";
for($j=1;$j<$num_subdirec;$j++)
{
	print WRBUFF $j;
	for($i=0;$i<@seg_edges;$i++)
	{
		$file1 = "$maindir/$seg_edges[$i]/1/$j/markers_segment$j.txt";
		$k = $j+1;
		$file2 = "$maindir/$seg_edges[$i]/1/$k/markers_segment$k.txt";
		$val1 = `tail -1 $file1|cut -f2 -d " "`;
		$val2 = `head -1 $file2|cut -f2 -d " "`;
		$val = ($val1+$val2)/2;	
		print  WRBUFF "\t$val";
	}
	print WRBUFF "\n";
}












##########################################################SUBROUTINES##########################################################################
sub median {
        @_ == 1 or die ('Sub usage: $median = median(\@array);');
        my ($array_ref) = @_;
        my $count = scalar @$array_ref;
        # Sort a COPY of the array, leaving the original untouched
        my @array = sort { $a <=> $b } @$array_ref;
        if ($count % 2)
        {
                return $array[int($count/2)];
        }
        else
        {
                return ($array[$count/2] + $array[$count/2 - 1]) / 2;
        }
}

sub mean {
@_ == 1 or die ('Sub usage: mean(\@array);');
my ($array_ref) = @_;
my $sum;
my $count = scalar @$array_ref;
foreach (@$array_ref) { $sum += $_; }
return $sum / $count;
}



sub sum {
@_ == 1 or die ('Sub usage: sum(\@array);');
my ($array_ref) = @_;
my $sum;
foreach (@$array_ref) { $sum += $_; }
return $sum ;
}

sub per1diff {
@_ == 1 or die ('Sub usage: per1diff(\@array);');
my ($array_ref) = @_;
my $sum=0;
my $count=0;
foreach (@$array_ref) {
        if($_ > 0.1)
        {
        $sum += $_;
        $count++;
        }
}
if($count == 0)
{
        return $sum;
}
else
{
        return $sum / $count;
}
}

sub per2diff {
@_ == 1 or die ('Sub usage: per2diff(\@array);');
my ($array_ref) = @_;
my $sum=0;
my $count=0;
foreach (@$array_ref) {
        if($_ > 0.2)
        {
        $sum += $_;
        $count++;
        }
}
if($count == 0)
{
        return $sum;
}
else
{
        return $sum / $count;
}
}

sub per5diff {
@_ == 1 or die ('Sub usage: per5diff(\@array);');
my ($array_ref) = @_;
my $sum=0;
my $count=0;
foreach (@$array_ref) {
        if($_ > 0.5)
        {
        $sum += $_;
        $count++;
        }
}
if($count == 0)
{
        return $sum;
}
else
{
        return $sum / $count;
}
}

